setwd('~/Research Papers/China CT-JME/china_ctjme_data')

library(foreign)
library(survival)
library(KMsurv)
library(MASS)
library(stargazer)
library(pastecs)
library(dplyr)
library(tidyr)
library(rmarkdown)
library(knitr)
library(ggplot2)
library(survminer)
library(pastecs)
library(pscl)
library(simPH)
library(ggthemes)
library(margins)
library(lubridate)
library(countrycode)
library(reshape)
library(stringr)
library(sandwich)
library(lmtest)
library(pscl)
library(readxl)
library(cowplot)
library(latex2exp)
library(readstata13)
library(Zelig)
library(car)
library(writexl)
library(readxl)
library(sjPlot)
library(sjmisc)
library(cowplot)
################################################################################

########## Plot of types of JMEs Figure 1 ######################################
d7<- read.csv("pla_jme_types.csv")

black.bold.22.text <- element_text(face = "bold", color = "black", size = 26)
g4<- ggplot(d7,aes(x=cat,y=count,fill=type))+
  geom_bar(stat="identity",position="dodge")+#scale_fill_grey()+
  scale_fill_manual(name="Measures of JMEs",
                    values=c('gray42','gray70'),
                    labels=c("Total JME days",
                             "Total JMEs"))+
  xlab("")+scale_x_discrete(labels=c("Combat", "Combat support", 
                                     "Competition","Counterterrorism",
                                     "MOOTW"))+
  ylab("Count")+#ggtitle("PLA Joint Military Exercises: 2002-2016")+
  geom_text(aes(label =label),vjust=0.6,hjust=-0.2,size = 6,position=position_dodge(0.8))+
  theme(axis.text.y =black.bold.22.text ,
        axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
        axis.text.x = element_text(size=22),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        text = element_text(size=26),legend.title = element_blank(),
        legend.background = element_rect(fill="white",size=1, linetype="solid", colour ="black"),
        legend.position=c(0.68, 0.88),legend.key.size=unit(1, "cm"),legend.key = element_rect(fill="white"),
        legend.text = element_text(colour="black", size = 18),panel.spacing=unit(c(0,0,2,0), "cm"))+
  coord_flip()
g4  
################################################################################

########## Time-series plot Figure 2############################################
##Plot time-series of counter-terrorism exercises (yearly)
d3<- read.csv("china_jmes_category.csv")
d5<- filter(d3, cat=="Counter-Terrorism")
d5$count<- 1
names(d5)
d6<- group_by(d5, Year, cat)
d6<- summarise(d6, count=sum(count), duration=sum(duration))

d6$Year<- as.numeric(d6$Year)



gts<- ggplot(d6,aes(x=Year,y=duration))+
  geom_line(colour="black", size=1.2)+#scale_fill_grey()+
  #scale_fill_manual(values=c('grey72','gray48','grey28'))+
  xlab("Year")+scale_x_continuous(breaks=seq(2002, 2016, 2))+
  ylab("Total number of CT-JME days")+#ggtitle("PLA Joint Military Exercises: 2002-2016")+
  theme(
    panel.background = element_blank(),
    axis.line = element_line(colour = "black"),
    panel.border = element_blank(),
    axis.title.y = element_text(margin = margin(t = 1, r = 10, b = 0, l = 0)),
    axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)),
    axis.text.x = element_text(size=18),
    axis.text.y = element_text(size=18),
    text = element_text(size=20),legend.title = element_blank(),
    legend.background = element_rect(fill="white",size=0.4, linetype="solid", colour ="black"),
    legend.position=c(0.1, 0.8),legend.key.size=unit(0.8, "cm"),legend.key = element_rect(fill="white"))
gts
################################################################################

### Codes for producing Model 1-6 in Table 1 ###################################
d2<- read.csv("CT_JME_v7.csv")
names(d2)

lg0 <- zeroinfl(ct_day ~ sqrtattack+as.factor(year)|
                  sqrtattack+as.factor(year),dist="negbin",link="logit",
                data=d2)
summary(lg0)
coeftest(lg0,vcov = vcovCL(lg0, cluster=d2$Country_Name))
rs0<- coeftest(lg0,vcov = vcovCL(lg0, cluster=d2$Country_Name))



lg1 <- zeroinfl(ct_day ~ lfdifull+as.factor(year)|
                  lfdifull+as.factor(year),dist="negbin",link="logit",
                data=d2)
summary(lg1)
coeftest(lg1,vcov = vcovCL(lg1, cluster=d2[which(d2$lfdifull!="NA"),]$Country_Name))
rs1<- coeftest(lg1,vcov = vcovCL(lg1, cluster=d2[which(d2$lfdifull!="NA"),]$Country_Name))


lg2 <- zeroinfl(ct_day ~ sqrtattack*lfdifull+as.factor(year)|
                  sqrtattack*lfdifull+as.factor(year),dist="negbin",link="logit",
                data=d2)
summary(lg2)
coeftest(lg2,vcov = vcovCL(lg2, cluster=d2[which(d2$lfdifull!="NA"),]$Country_Name))
rs2<- coeftest(lg2,vcov = vcovCL(lg2, cluster=d2[which(d2$lfdifull!="NA"),]$Country_Name))


lg4 <- zeroinfl(ct_day ~ sqrtattack*lfdifull+
                  log(gdp_wb)+larmtrade+polity2+us_defense+sco+extradition+
                  log(distcap)+obor+as.factor(year)|
                  sqrtattack*lfdifull+
                  log(gdp_wb)+larmtrade+polity2+us_defense+sco+extradition+
                  log(distcap)+obor+as.factor(year),dist="negbin",link="logit",
                data=d2)
summary(lg4)
coeftest(lg4,vcov = vcovCL(lg4, cluster=d2[which(d2$lfdifull!="NA" & d2$polity2!="NA"),]$Country_Name))
rs4<- coeftest(lg4,vcov = vcovCL(lg4, cluster=d2[which(d2$lfdifull!="NA" & d2$polity2!="NA"),]$Country_Name))


d2$sqrtdom<- sqrt(d2$domestic)
d2$sqrtint<- sqrt(d2$int_any)
lg5 <- zeroinfl(ct_day ~ sqrtdom*lfdifull+sqrtint*lfdifull+
                  as.factor(year)|
                  sqrtdom*lfdifull+sqrtint*lfdifull+
                  as.factor(year),dist="negbin",link="logit",
                data=d2)
summary(lg5)
coeftest(lg5,vcov = vcovCL(lg5, cluster=d2[which(d2$lfdifull!="NA"),]$Country_Name))
rs5<- coeftest(lg5,vcov = vcovCL(lg5, cluster=d2[which(d2$lfdifull!="NA"),]$Country_Name))

lg6 <- zeroinfl(ct_day ~ sqrtdom*lfdifull+sqrtint*lfdifull+
                  log(gdp_wb)+larmtrade+polity2+us_defense+sco+extradition+
                  log(distcap)+obor+as.factor(year)|
                  sqrtdom*lfdifull+sqrtint*lfdifull+
                  log(gdp_wb)+larmtrade+polity2+us_defense+sco+extradition+
                  log(distcap)+obor+as.factor(year),dist="negbin",link="logit",
                data=d2)
summary(lg6)
coeftest(lg6,vcov = vcovCL(lg6, cluster=d2[which(d2$lfdifull!="NA" & d2$polity2!="NA"),]$Country_Name))
rs6<- coeftest(lg6,vcov = vcovCL(lg6, cluster=d2[which(d2$lfdifull!="NA" & d2$polity2!="NA"),]$Country_Name))


stargazer(rs0, rs1,rs2,rs4,rs5,rs6,
          align=TRUE, no.space = TRUE,
          column.sep.width = "-16pt")
################################################################################

## Codes for producing Figure 4#################################################
lg0 <- zeroinfl(ct_day ~ sqrtattack*lfdifull+
                  log(gdp_wb)+larmtrade+polity2+us_defense+sco+extradition+
                  log(distcap)+obor+as.factor(year)|
                  sqrtattack*lfdifull+
                  log(gdp_wb)+larmtrade+polity2+us_defense+sco+extradition+
                  log(distcap)+obor+as.factor(year),dist="negbin",link="logit",
                data=d2)
summary(lg0)

tenvalues <- seq(0, 62, by=3)
count.0 <- sapply(tenvalues, FUN=function(x){
  mean(predict(lg0, type = "response", 
               newdata = mutate(d2,
                                sqrtattack=x, 
                                lfdifull=mean(d2$lfdifull,na.rm=TRUE)+
                                  2*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})
count.1 <- sapply(tenvalues, FUN=function(x){
  mean(predict(lg0, type = "response", 
               newdata = mutate(d2,
                                sqrtattack=x, 
                                lfdifull=mean(d2$lfdifull,na.rm=TRUE)+
                                  1*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})

count.2 <- sapply(tenvalues, FUN=function(x){
  mean(predict(lg0, type = "response", 
               newdata = mutate(d2,
                                sqrtattack=x, 
                                lfdifull=mean(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})

count.3 <- sapply(tenvalues, FUN=function(x){
  mean(predict(lg0, type = "response", 
               newdata = mutate(d2,
                                sqrtattack=x, 
                                lfdifull=mean(d2$lfdifull,na.rm=TRUE)-
                                  1*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})

QI <- data.frame(sqrtattack = tenvalues,
                 predcount = c(count.0, count.1, count.2, count.3),
                 directchange = rep(c("High FDI",
                                      "Moderate FDI",
                                      "Mean FDI", 
                                      "Low FDI"), each=21)
)

B <- coef(lg0)
V <- vcovCL(lg0,vcov = vcovCL(lg0, cluster=d2[which(d2$lfdifull!="NA" & d2$polity2!="NA"),]$Country_Name))

sim.coefs <- mvrnorm(1000, mu=B, Sigma=V)


zinb2 <- lg0
zinb2$coefficients<- unlist(zinb2$coefficients) #since the coefficients of a ZINB model has two parts, we need first to make it a single matrix by unlisting it
sim.qi <- apply(sim.coefs, 1, FUN=function(y){
  zinb2$coefficients<- y
  zinb2$coefficients<- list(count=zinb2$coefficients[1:25],
                            zero=zinb2$coefficients[26:50]) #re-list the two set of coefficients so that R could recgonize it as a ZINB object
  count.0 <- sapply(tenvalues, FUN=function(x){
    mean(predict(zinb2, type = "response", 
                 newdata = mutate(d2,
                                  sqrtattack=x, 
                                  lfdifull=mean(d2$lfdifull,na.rm=TRUE)+
                                    2*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})
  
  count.1 <- sapply(tenvalues, FUN=function(x){
    mean(predict(zinb2, type = "response", 
                 newdata = mutate(d2,
                                  sqrtattack=x, 
                                  lfdifull=mean(d2$lfdifull,na.rm=TRUE)+
                                    1*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})
  
  count.2 <- sapply(tenvalues, FUN=function(x){
    mean(predict(zinb2, type = "response", 
                 newdata = mutate(d2,
                                  sqrtattack=x, 
                                  lfdifull=mean(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})
  
  count.3 <- sapply(tenvalues, FUN=function(x){
    mean(predict(zinb2, type = "response", 
                 newdata = mutate(d2,
                                  sqrtattack=x, 
                                  lfdifull=mean(d2$lfdifull,na.rm=TRUE)-
                                    1*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})
  
  return(c(count.0, count.1, count.2, count.3))  
})
sim.qi <- t(sim.qi)

predcount.se <- apply(sim.qi, 2, sd)
predcount.lb <- apply(sim.qi, 2, FUN=function(x){
  quantile(x, .025)
})
predcount.ub <- apply(sim.qi, 2, FUN=function(x){
  quantile(x, .975)
})
QI2 <- mutate(QI, 
              predcount.se = predcount.se,
              predcount.lb = predcount.lb,
              predcount.ub = predcount.ub)

QI2$attack<- QI2$sqrtattack*QI2$sqrtattack

QI2$types<- rep(c("Group 1",
                  "Group 2", 
                  "Group 3",
                  "Group 4"), each=21)

g1 <- ggplot(QI2[which(QI2$attack<=400),],
             aes(x = sqrtattack, y = predcount, group=types))+
  #facet_grid(. ~ directchange)+
  geom_line(aes(linetype=types, size=types))+
  scale_size_manual(values=c(1, 1, 1, 1), labels=c("High FDI (+2 sd)", "Moderately High FDI (+1 sd)", "Mean FDI (mean)", "Low FDI (-1 sd)"))+
  scale_linetype_manual(values=c("solid","dotted","dashed", "dotdash"), labels=c("High FDI (+2 sd)", "Moderately High FDI (+1 sd)", "Mean FDI (mean)", "Low FDI (-1 sd)"))+
  geom_ribbon(aes(ymin=predcount.lb,ymax=predcount.ub),alpha=0.2) +
  geom_hline(yintercept=0, linetype="solid", color = "red",size=1)+
  ylab("Number of CT-JME Days") +
  xlab("Terrorist Attack (square root value)") +
  #scale_x_continuous(breaks = round(seq(0, 16, by = 2),1))+
  theme_bw()+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 28, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 30, r = 0, b = 0, l = 0)))+
  theme(text = element_text(size=30),legend.title = element_blank(),
        strip.text.x = element_text(size = 30),
        legend.background = element_rect(fill="white",size=0.6, linetype="solid", colour ="black"),
        legend.text=element_text(size=30),
        legend.position=c(0.2, 0.8),legend.key.size=unit(1.5, "cm"),legend.key = element_rect(fill="white"),
        plot.caption = element_text(size = 25, hjust = 0,vjust = -3, face="italic"),
        plot.margin = margin(1, 2, 1, 2, "cm"))+
  labs(caption = "Note: shaded areas represent 95% credible intervals generated through 1000 draws of new coefficients from the posterior.")
g1
################################################################################

## Codes for producing Figure 5#################################################
lg0 <- zeroinfl(ct_day ~ sqrtattack*lfdifull+
                  log(gdp_wb)+larmtrade+polity2+us_defense+sco+extradition+
                  log(distcap)+obor+as.factor(year)|
                  sqrtattack*lfdifull+
                  log(gdp_wb)+larmtrade+polity2+us_defense+sco+extradition+
                  log(distcap)+obor+as.factor(year),dist="negbin",link="logit",
                data=d2)
summary(lg0)

tenvalues <- seq(4.6, 20, by=0.6)

B <- coef(lg0)
V <- vcovCL(lg0,vcov = vcovCL(lg0, cluster=d2[which(d2$lfdifull!="NA" & d2$polity2!="NA"),]$Country_Name))
sim.coefs <- mvrnorm(1000, mu=B, Sigma=V)

zinb1 <- lg0
zinb1$coefficients<- unlist(zinb1$coefficients)
QI <- sapply(tenvalues, FUN=function(pval){
  count.pred <- predict(lg0, type="response", newdata=mutate(d2, lfdifull=pval))
  prob.zero <- predict(lg0, type="zero", newdata=mutate(d2, lfdifull=pval))
  y.star <- log(count.pred/(lg0$theta*(1-prob.zero)))
  coef.attack <- coef(lg0)[["count_sqrtattack"]]
  coef.inter<- coef(lg0)[["count_sqrtattack:lfdifull"]]
  marg.prob <- lg0$theta *(coef.attack+pval*coef.inter)*exp(y.star)
  point.est <- mean(marg.prob, na.rm=TRUE)
  marg.prob.sim <- apply(sim.coefs, 1, FUN=function(x){
    zinb1$coefficients <- x
    zinb1$coefficients <- list(count=zinb1$coefficients[1:25], zero=zinb1$coefficients[26:50])
    count.pred <- predict(zinb1, type="response", newdata=mutate(d2, lfdifull=pval))
    prob.zero <- predict(zinb1, type="zero", newdata=mutate(d2, lfdifull=pval))
    y.star <- log(count.pred/(zinb1$theta*(1-prob.zero)))
    coef.attack <- zinb1$coefficients$count[[2]]
    coef.inter <- zinb1$coefficients$count[[25]]
    marg.prob <- zinb1$theta*(coef.attack+pval*coef.inter)*exp(y.star)
    return(mean(marg.prob, na.rm=TRUE))
  })
  se <- sd(marg.prob.sim)
  lb <- quantile(marg.prob.sim, .025)
  ub <- quantile(marg.prob.sim, .975)
  return(c(point.est, se, lb, ub))
})
QI <- data.frame(tenvalues, t(QI))
colnames(QI) <- c("lfdi", "point.est", "se", "lb", "ub")

QI$fdi<- exp(QI$lfdi)

g2 <- ggplot(QI,aes(x = lfdi, y = point.est))+
  #facet_grid(. ~ directchange)+
  geom_line()+
  scale_size_manual(values=c(1))+
  scale_linetype_manual(values=c("solid"))+
  geom_ribbon(aes(ymin=lb,ymax=ub),alpha=0.2) +
  geom_hline(yintercept=0, linetype="solid", color = "red",size=1)+
  ylab("Marginal Effect of Terrorist Attack on CT-JME Days") +
  xlab("FDI (log value)") +
  #scale_x_continuous(breaks = round(seq(0, 16, by = 2),1))+
  theme_bw()+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 28, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 30, r = 0, b = 0, l = 0)))+
  theme(text = element_text(size=30),legend.title = element_blank(),
        strip.text.x = element_text(size = 30),
        legend.background = element_rect(fill="white",size=0.5, linetype="solid", colour ="black"),
        legend.text=element_text(size=30),
        legend.position=c(0.1, 0.8),legend.key.size=unit(0.8, "cm"),legend.key = element_rect(fill="white"),
        plot.caption = element_text(size = 25, hjust = 0,vjust = -3, face="italic"),
        plot.margin = margin(1, 2, 1, 2, "cm"))+
  labs(caption = "Note: shaded areas represent 95% credible intervals generated through 1000 draws of new coefficients from the posterior.")
g2
################################################################################

## Codes for producing Figure 6#################################################
lg0 <- zeroinfl(ct_day ~ sqrtdom*lfdifull+sqrtint*lfdifull+
                  log(gdp_wb)+larmtrade+polity2+us_defense+sco+extradition+
                  log(distcap)+obor+as.factor(year)|
                  sqrtdom*lfdifull+sqrtint*lfdifull+
                  log(gdp_wb)+larmtrade+polity2+us_defense+sco+extradition+
                  log(distcap)+obor+as.factor(year),dist="negbin",link="logit",
                data=d2)
summary(lg0)
coeftest(lg0,vcov = vcovCL(lg0, cluster=d2[which(d2$lfdifull!="NA" & d2$polity2!="NA"),]$Country_Name))

## Domestic
tenvalues <- seq(0, 34, by=1)
count.0 <- sapply(tenvalues, FUN=function(x){
  mean(predict(lg0, type = "response", 
               newdata = mutate(d2,
                                sqrtdom=x, 
                                sqrtint=mean(d2$sqrtint),
                                lfdifull=mean(d2$lfdifull,na.rm=TRUE)+
                                  2*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})

count.1 <- sapply(tenvalues, FUN=function(x){
  mean(predict(lg0, type = "response", 
               newdata = mutate(d2,
                                sqrtdom=x, 
                                sqrtint=mean(d2$sqrtint),
                                lfdifull=mean(d2$lfdifull,na.rm=TRUE)+
                                  1*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})

count.2 <- sapply(tenvalues, FUN=function(x){
  mean(predict(lg0, type = "response", 
               newdata = mutate(d2,
                                sqrtdom=x, 
                                sqrtint=mean(d2$sqrtint),
                                lfdifull=mean(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})

count.3 <- sapply(tenvalues, FUN=function(x){
  mean(predict(lg0, type = "response", 
               newdata = mutate(d2,
                                sqrtdom=x, 
                                sqrtint=mean(d2$sqrtint),
                                lfdifull=mean(d2$lfdifull,na.rm=TRUE)-
                                  1*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})

QI <- data.frame(sqrtattack = tenvalues,
                 predcount = c(count.0, count.1, count.2, count.3),
                 directchange = rep(c("High FDI",
                                      "Moderate FDI",
                                      "Mean FDI", 
                                      "Low FDI"), each=35)
)

B <- coef(lg0)
V <- vcovCL(lg0,vcov = vcovCL(lg0, cluster=d2[which(d2$lfdifull!="NA" & d2$polity2!="NA"),]$Country_Name))

sim.coefs <- mvrnorm(1000, mu=B, Sigma=V)


zinb2 <- lg0
zinb2$coefficients<- unlist(zinb2$coefficients) #since the coefficients of a ZINB model has two parts, we need first to make it a single matrix by unlisting it
sim.qi <- apply(sim.coefs, 1, FUN=function(y){
  zinb2$coefficients<- y
  zinb2$coefficients<- list(count=zinb2$coefficients[1:27],
                            zero=zinb2$coefficients[28:54]) #re-list the two set of coefficients so that R could recgonize it as a ZINB object
  count.0 <- sapply(tenvalues, FUN=function(x){
    mean(predict(zinb2, type = "response", 
                 newdata = mutate(d2,
                                  sqrtdom=x, 
                                  sqrtint=mean(d2$sqrtint),
                                  lfdifull=mean(d2$lfdifull,na.rm=TRUE)+
                                    2*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})
  
  count.1 <- sapply(tenvalues, FUN=function(x){
    mean(predict(zinb2, type = "response", 
                 newdata = mutate(d2,
                                  sqrtdom=x, 
                                  sqrtint=mean(d2$sqrtint),
                                  lfdifull=mean(d2$lfdifull,na.rm=TRUE)+
                                    1*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})
  
  count.2 <- sapply(tenvalues, FUN=function(x){
    mean(predict(zinb2, type = "response", 
                 newdata = mutate(d2,
                                  sqrtdom=x, 
                                  sqrtint=mean(d2$sqrtint),
                                  lfdifull=mean(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})
  
  count.3 <- sapply(tenvalues, FUN=function(x){
    mean(predict(zinb2, type = "response", 
                 newdata = mutate(d2,
                                  sqrtdom=x, 
                                  sqrtint=mean(d2$sqrtint),
                                  lfdifull=mean(d2$lfdifull,na.rm=TRUE)-
                                    1*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})
  
  return(c(count.0, count.1, count.2, count.3))  
})
sim.qi <- t(sim.qi)

predcount.se <- apply(sim.qi, 2, sd)
predcount.lb <- apply(sim.qi, 2, FUN=function(x){
  quantile(x, .025)
})
predcount.ub <- apply(sim.qi, 2, FUN=function(x){
  quantile(x, .975)
})
QI3 <- mutate(QI, 
              predcount.se = predcount.se,
              predcount.lb = predcount.lb,
              predcount.ub = predcount.ub)

QI3$attack<- QI3$sqrtattack*QI3$sqrtattack

QI3$types<- rep(c("Group 1",
                  "Group 2", 
                  "Group 3",
                  "Group 4"), each=35)
QI3$category<- "Domestic Attack"

## Transnational
tenvalues <- seq(0, 34, by=1)
count.0 <- sapply(tenvalues, FUN=function(x){
  mean(predict(lg0, type = "response", 
               newdata = mutate(d2,
                                sqrtdom=mean(d2$sqrtdom),
                                sqrtint=x, 
                                lfdifull=mean(d2$lfdifull,na.rm=TRUE)+
                                  2*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})

count.1 <- sapply(tenvalues, FUN=function(x){
  mean(predict(lg0, type = "response", 
               newdata = mutate(d2,
                                sqrtdom=mean(d2$sqrtdom),
                                sqrtint=x, 
                                lfdifull=mean(d2$lfdifull,na.rm=TRUE)+
                                  1*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})


count.2 <- sapply(tenvalues, FUN=function(x){
  mean(predict(lg0, type = "response", 
               newdata = mutate(d2,
                                sqrtdom=mean(d2$sqrtdom),
                                sqrtint=x, 
                                lfdifull=mean(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})

count.3 <- sapply(tenvalues, FUN=function(x){
  mean(predict(lg0, type = "response", 
               newdata = mutate(d2,
                                sqrtdom=mean(d2$sqrtdom),
                                sqrtint=x, 
                                lfdifull=mean(d2$lfdifull,na.rm=TRUE)-
                                  1*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})

QI <- data.frame(sqrtattack = tenvalues,
                 predcount = c(count.0, count.1, count.2, count.3),
                 directchange = rep(c("High FDI",
                                      "Moderate FDI",
                                      "Mean FDI", 
                                      "Low FDI"), each=35)
)

B <- coef(lg0)
V <- vcovCL(lg0,vcov = vcovCL(lg0, cluster=d2[which(d2$lfdifull!="NA" & d2$polity2!="NA"),]$Country_Name))

sim.coefs <- mvrnorm(1000, mu=B, Sigma=V)


zinb2 <- lg0
zinb2$coefficients<- unlist(zinb2$coefficients) #since the coefficients of a ZINB model has two parts, we need first to make it a single matrix by unlisting it
sim.qi <- apply(sim.coefs, 1, FUN=function(y){
  zinb2$coefficients<- y
  zinb2$coefficients<- list(count=zinb2$coefficients[1:27],
                            zero=zinb2$coefficients[28:54]) #re-list the two set of coefficients so that R could recgonize it as a ZINB object
  count.0 <- sapply(tenvalues, FUN=function(x){
    mean(predict(zinb2, type = "response", 
                 newdata = mutate(d2,
                                  sqrtdom=mean(d2$sqrtdom),
                                  sqrtint=x, 
                                  lfdifull=mean(d2$lfdifull,na.rm=TRUE)+
                                    2*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})
  
  count.1 <- sapply(tenvalues, FUN=function(x){
    mean(predict(zinb2, type = "response", 
                 newdata = mutate(d2,
                                  sqrtdom=mean(d2$sqrtdom),
                                  sqrtint=x, 
                                  lfdifull=mean(d2$lfdifull,na.rm=TRUE)+
                                    1*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})
  
  
  count.2 <- sapply(tenvalues, FUN=function(x){
    mean(predict(zinb2, type = "response", 
                 newdata = mutate(d2,
                                  sqrtdom=mean(d2$sqrtdom),
                                  sqrtint=x, 
                                  lfdifull=mean(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})
  
  count.3 <- sapply(tenvalues, FUN=function(x){
    mean(predict(zinb2, type = "response", 
                 newdata = mutate(d2,
                                  sqrtdom=mean(d2$sqrtdom),
                                  sqrtint=x, 
                                  lfdifull=mean(d2$lfdifull,na.rm=TRUE)-
                                    1*sd(d2$lfdifull,na.rm=TRUE))), na.rm=TRUE)})
  
  return(c(count.0, count.1, count.2, count.3))  
})
sim.qi <- t(sim.qi)

predcount.se <- apply(sim.qi, 2, sd)
predcount.lb <- apply(sim.qi, 2, FUN=function(x){
  quantile(x, .025)
})
predcount.ub <- apply(sim.qi, 2, FUN=function(x){
  quantile(x, .975)
})
QI4 <- mutate(QI, 
              predcount.se = predcount.se,
              predcount.lb = predcount.lb,
              predcount.ub = predcount.ub)

QI4$attack<- QI4$sqrtattack*QI3$sqrtattack
QI4$types<- rep(c("Group 1",
                  "Group 2", 
                  "Group 3",
                  "Group 4"), each=35)
QI4$category<- "Transnational Attack"

QI5<- rbind(QI3, QI4)



g1 <- ggplot(QI5[which(QI5$attack<=100),],
             aes(x = sqrtattack, y = predcount, group=types))+facet_grid(. ~ category)+
  #facet_grid(. ~ directchange)+
  geom_line(aes(linetype=types, size=types))+
  scale_size_manual(values=c(1, 1, 1, 1), labels=c("High FDI (+2 sd)", "Moderately High FDI (+1 sd)", "Mean FDI (mean)", "Low FDI (-1 sd)"))+
  scale_linetype_manual(values=c("solid","dotted","dashed", "dotdash"), labels=c("High FDI (+2 sd)", "Moderately High FDI (+1 sd)", "Mean FDI (mean)", "Low FDI (-1 sd)"))+
  geom_ribbon(aes(ymin=predcount.lb,ymax=predcount.ub),alpha=0.2) +
  geom_hline(yintercept=0, linetype="solid", color = "red",size=1)+
  ylab("Number of CT-JME Days") +
  xlab("Terrorist Attack (square root value)") +
  #scale_x_continuous(breaks = round(seq(0, 16, by = 2),1))+
  theme_bw()+
  theme(axis.title.y = element_text(margin = margin(t = 30, r = 28, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 30, r = 0, b = 0, l = 0)))+
  theme(text = element_text(size=30),legend.title = element_blank(),
        strip.text.x = element_text(size = 30),
        legend.background = element_rect(fill="white",size=0.6, linetype="solid", colour ="black"),
        legend.text=element_text(size=30),
        legend.position="right",legend.key.size=unit(1.5, "cm"),legend.key = element_rect(fill="white"),
        plot.caption = element_text(size = 25, hjust = 0,vjust = -3, face="italic"),
        plot.margin = margin(1, 2, 1, 2, "cm"))+
  labs(caption = "Note: shaded areas represent 95% credible intervals generated through 1000 draws of new coefficients from the posterior.")
g1
################################################################################